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ABSTRACT 

The HI distribution at the epoch of reionization (EOR) is largely determined by the 
sizes and distribution of the ionized regions. In the scenario where the ionized regions 
\ have comoving sizes of the order of a few Mpc, the large scale statistical properties of 

^ ' the HI distribution are dominated by the Poisson noise of the discrete ionized regions, 

and it is highly non-Gaussian. We investigate the possibility of probing reionization 
by studying these non-Gaussian features using future radio interferometric observa- 
tions of redshifted 21 cm HI radiation. We develop a formalism relating correlations 
between the visibilities measured at three different baselines and frequencies to the 
\ bispectrum of HI fluctuations. For visibilities at the same frequency, this signal is 

■ found to be of the same order as the two visibility correlation which probes the HI 

^ (— | power spectrum. For visibilities at different frequencies, we find that the correlations 

O decay within a frequency difference of ~ 1 MHz. This implies that it is, in principle, 

straightforward to extract this HI signal from various contaminants which are believed 
£5 ' to have a continuum spectra and are expected to be correlated even at large frequency 

^ \ separations. 



Key words: cosmology: theory - cosmology: large scale structure of universe - inter- 
galactic medium - diffuse radiation 



1 INTRODUCTION: 

There has recently been a lot of interest in understanding exactly how and when the universe was reionized. There now are 
significant observational constraints mainly from three different kinds of observations. The observation of quasars at redshift 
z ~ 6 which show strong HI absorption (Becker et al. 2001) indicates that at least 1% of the total hydrogen mass at z ~ 6 is 
neutral (Fan et al. 2002), and the neutral mass fraction decreases rapidly at lower redshifts. This is a strong indication that 
the epoch of reionization ended at z ~ 6 . Observations of the CMBR polarization, generated through Thomson scattering of 
CMBR photons by free electrons along the line of sight, indicates that the reionization began at a redshift z > 14. On the 
other hand, the observed anisotropies of the CMBR indicate that the total optical depth of the Thomson scattering is not 
extremely high, suggesting that reionization could not have started at redshift much higher than about 30 (Kogut et al. 2003; 
Spergel et al. 2003). A third constraint comes from determinations of the IGM temperature from observations of the Lya 
forest in the z range 2 to 4 which indicates a complex reionization history with there possibly being an order unity change in 
the neutral hydrogen fraction at z < 10 (Theuns et al. 2002; Hui & Haiman 2003). 

Mapping the HI distribution at high redshifts using radio observations of the redshifted 21 cm radiation (Madau, Meiksin 
& Rees 1997; Scott & Rees 1990; Kumar, Padmanabhan & Subramanian 1985 ) holds the possibility of probing the transition 
from a largely neutral to a largely ionized universe at a level of detail surpassing any other techniques. Zaldarriaga, Furlanetto 
& Hernquist (2003) (hereafter ZFH) have developed a statistical technique based on the angular power spectrum, on lines 
similar to the analysis of CMBR anisotropies, for analysing the HI signal from the epoch of reionization (EOR) in radio 
interferometric observations. Extracting the HI signal from various Galactic and extra Galactic contaminants (eg. Cooray & 
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Furlanotto 2004; DiMatteo. et al. 2004; Gnedin & Shaver 2003; Oh & Mack 2003; DiMatteo. et al. 2002; Shaver et al. 1999; ) 
is one of the most important challenges. Most of the known contaminants are expected to have continuum spectra, and ZFH 
show that it should in principle be possible to extract the HI signal using the fact that, unlike the contaminants, it will be 
uncorrelated at two slightly different frequencies. The frequency dependence of the angular power spectrum of the HI signal 
and foregrounds has recently been analysed in detail by Santos, Cooray & Knox (2004) . 

An alternative statistical technique for analysing the HI signal is to study the correlations between the complex visibilities 
measured at different baselines and frequencies in radio-interferometric observations. This has been developed in the context 
of observing HI from z < 6 (Bharadwaj & Sethi 2001; Bharadwaj & Pandey 2003; Bharadwaj & Srikant 2004) and later 
generalized to the EOR signal in Bharadwaj & Ali 2004 (hereafter BA). The possibility of using visibility correlations to 
quantify the EOR signal has also been proposed by Morales & Hewitt (2003) who further discuss how the different frequency 
signatures of the contaminants and the HI signal can be used to distinguish between the two. Recently (Morales 2004) has 
addressed the issue of the power spectrum sensitivity of the EOR HI signal. 

Various investigations (eg. ZFH, BA) show that the power spectrum of HI fluctuations at EOR has contributions from 
mainly two distinct effects, the clustering of the hydrogen which, on large scales, is assumed to follow the dark matter 
distribution and the fluctuations arising from the presence of discrete regions of ionized gas surrounding the sources responsible 
for reionizing the universe. The details of the reionization process are not very well understood (eg. Barkana & Loeb 2001), 
and the shape, size and distribution of these ionized regions is one of the very important issues which will be probed by 21 
cm HI observations. There has recently been progress in analytically modeling the growth of the ionized regions (Furlanetto, 
Zaldarriaga & Hernquist 2004a) (hereafter FZH) based on the findings of simulations (Ciardi et al. 2003; Sokasian et al. 2003a; 
Sokasian et al. 2003b; Nusser et al. 2002; Benson et al. 2001; Gnedin 2000) which show that there will not be a large number 
of small HII regions around individual ionizing sources, rather there will be a few large ionized regions centered on places 
where the ionized sources are clustered. The size of these ionized regions are expected to be around a few Mpc (comoving) or 
possibly larger at EOR. In such a scenarios, on scales larger than the size of the individual ionized regions, the HI signal will 
be dominated by the Poisson noise arising from the discrete nature of the ionized regions (eg. ZFH, BA, FZH, Furlanetto, 
Zaldarriaga & Hernquist 2004b). Further, the HI signal is expected to be highly non-Gaussian . 

Nearly all of the work on quantifying the EOR HI signal expected in radio interferometric observations has focused 
on the two point statistics namely the angular power spectrum and the correlations between pairs of visibilities. Both these 
quantities are actually equivalent and they basically probe the power spectrum of HI fluctuations at EOR. The power spectrum 
completely quantifies a Gaussian random field, but the higher order statistics would contain independent information if the HI 
fluctuations at EOR were not a Gaussian random field. FZH have used the pixel distribution function, a one-point statistics, 
to quantify non-Gaussian features in the HI distribution. He et al. (2004) have studied the non-Gaussian features that arise 
in the HI distribution in the log-normal model. 

In this paper we address the issue of quantifying the non-Gaussian features of the HI signal expected in radio interfer- 
ometric observations. In particular, we focus on the correlation between three visibilities. This is expected to be zero if the 
signal were a Gaussian random field, and deviations from zero are a clear signature of the non-Gaussian properties of the HI 
distribution. Here we derive the relation between the three visibility correlation and the bispectrum of the HI fluctuations. 
The bispectrum quantifies correlations between three Fourier modes, and this is non-zero only when there are phase corre- 
lations between different modes. The three visibility correlation, as we show, is comparable to the correlations between two 
visibilities and this leads us to speculate that this will play an important role in detecting the HI signal. Further, the higher 
order correlations contains independent information, and observing these would throw independent light on the topology and 
morphology of the HI distribution at EOR. 

Finally, an outline of the paper. In Section 2. we present the formalism relating the three visibility correlation to the HI 
bispectrum. In Section 3. we introduce a simple model for the HI distribution at reionization and calculate its bispectrum. In 
Section 4. we present results for the three visibility correlation expected from HI at reionization and discuss some consequences. 

It may also be noted that we use the values (ft m o, ^ao, fi&ft 2 , h) — (0.3,0.7,0.02,0.7) for the cosmological parameters 
thoughout. 



2 FORMALISM FOR THREE VISIBILITY CORRELATION 

In this section we follow the notation used in BA which also contains a more detailed discussion of the formalism for 
calculating the HI signal. The HI radiation at frequency 1420 MHz in the rest frame of the hydrogen is redshifted to a 
frequency v = 1420/(1 + z) MHz for an observer at present. The expansion of the universe and the HI peculiar velocity both 
contribute to the redshift. Incorporating these effects, the specific intensity I v (n) of redshifted 21 cm HI radiation at frequency 
v and direction n can be written as I v (n) = I v {z) x ??Hi(n, z) where 




(1) 
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and 



x Phi f, T 7 \ (1 + 0) dv . . 

It should be noted that the terms on the right hand side of equations (1) and (2) refer to the epoch when the HI radiation 
originated. Here H(z) the Hubble parameter, pn the mean cosmological density of hydrogen and r (or r„) the comoving 
distance to the HI calculated ignoring peculiar velocities, depend only on z. The quantities pm the HI density, T 7 the CMBR 
temperature, T s the HI spin temperature and v the radial component of the HI peculiar velocity also vary with position and 
should be evaluated at x = r„n ie. the position where the radiation originated. It may be noted that ??hi(x, 2), the 21cm 
radiation efficiency, was originally introduced by Madau, Meiksin & Rees (1997) who did not include peculiar velocities. As 
shown in BA, equation (2) includes an extra term which arises when the effect of the HI peculiar velocities are included. The 
quantity r)m(n, z) incorporate the details of the HI evolution including effects of heating, reionization and density fluctuations 
due to structure formation. 

We next introduce r/m(k, z), the Fourier transform of J)m(y, 2), 

^m(y,2) = y ^-^e~' k y ?7Hi(k,2) . (3) 

where y refers to an arbitrary comoving position. Using this we can express »?Hi(n, 2) as 

Uhi(A,*) = J j^e-'^rhnfrz) (4) 

where it is understood that this refers to the position x = r v n. 

The ensemble average of various products of 77m (ki, 2) are used to quantify the statistical properties of the fluctuation 
in the HI distribution. We first consider the HI power spectrum PHi(ki, 2) defined through 

(r?Hi(ki,2)j7Hi(k 2 ,2)> = (2TvfSl(k 1 + k 2 ) P m (ki, 2) (5) 

where 8% is the three dimensional Dirac delta function. The power spectrum completely quantifies all properties of the HI 
distribution if the fluctuations are a Gaussian random field. The higher order statistics contain independent information if 
the fluctuations are not a Gaussian random field. Here we proceed one step beyond the power spectrum and also consider the 
HI bispectrum Pm(ki, k 2 , k3, 2) defined through 

(^Hi(ki,2)r7Hi(k 2 ,2)?7Hi(k3,2)) = (27r) 3 5|,(ki +k 2 + k 3 )BHi(ki,k 2 ,k 3 ,2). (6) 

We next mention a few well known properties of the power spectrum and bispectrum which are relevant to the discussion. 
The fact that not all modes are correlated, reflected in the Dirac delta functions in eq. (5) and (6), is a consequence of the 
assumption that HI fluctuations are statistical homogeneous. Further, Pm(k) is isotropic ie. does not depend on the direction 
of k, if the effects of the peculiar velocity are ignored. The redshift space distortion caused by the peculiar velocities breaks 
the isotropy of Pm(k) which now depends on the orientation of k with respect to the line of sight. Similarly, ignoring redshift 
space distortions, i?Hi(ki, k 2 , k.3) depends only on the triangle formed by the wave vectors ki, k 2 and k3, and this is completely 
specified by the magnitude of the three vectors (fei, fe 2 , £3). The bispectrum also depends on how the triangle is oriented with 
respect to the line of sight if redshift space distortions are included. Finally, we note that both the power-spectrum and the 
bispectrum are real quantities. While the power spectrum is necessarily positive, there is no such restriction on the bispectrum. 

We now shift our attention to radio interferometric observations of redshifted HI using an array of low frequency radio 
antennas distributed on a plane. The antennas all point in the same direction m which we take to be vertically up wards. The 
beam pattern A(6) quantifies how the individual antenna, pointing up wards, responds to signals from different directions in 
the sky. This is assumed to be a Gaussian A(ff) = e~ with 9o <§C 1 i.e. the beam width of the antennas (in radians) is 
small, and the part of the sky which contributes to the signal can be well approximated by a plane. In this approximation 
the unit vector n can be represented by n = m + 9, where 6 is a two dimensional vector in the plane of the sky. Using this 
the angular fluctuations in the specific intensity <5/„ can be expressed as 

«„(&) = /„/ ^ e -^(<=ii+^-^ H i(k,2) (7) 

where fcy = k • m and kx are respectively the components of k parallel and perpendicular to m. The component kj_ lies in 
the plane of the sky. 

The quantity measured in interferometric observations is the complex visibility V(XJ, v) which is recorded for every 
independent pair of antennas at every frequency channel in the band of observations. For any pair of antennas, U = d/A 
quantifies the separation d in units of the wavelength A, we refer to this dimensionless quantity U as a baseline. A typical 
radio interferometric array simultaneously measures visibilities at a large number of baselines and frequency channels, and 
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V(TJ,v)= I d 2 8A{8)I v {e)e~ t2lTV - e . (8) 



The visibilities record only the angular fluctuations in and the visibilities arising from angular fluctuations in the 

HI radiation are 

-/3;„ 

»fc||r„ (g) 



V{\J,v) = L J ^a(U-gk x )^i(k, Z )e- 



where a(!7) the Fourier transform of the antenna beam pattern A(6), which for a Gaussian beam A{9) = e B l e a gives 
the Fourier transform also to be a Gaussian a(U) = nO 2 , exp [-TT 2 6lU 2 ] which we use in the rest of this paper. 

In this paper we quantify the statistical properties of the quantity measured in radio-interferometric observations, namely 
the visibilities at different baselines and frequencies. Further, we study their relation to the statistical properties of the HI 
distribution. To this end, we introduce the notation 

S 2 (Ui, U 2 , Av) = {V(Vi,v + Av)V{U 2 , v)) (10) 
and 

S 3 (Ui, U„, U 3 , Avi, A1/2) = <V(Ui, v + Avi)V(U 2 , v + Av 2 ) V{U 3 , v)) (11) 

to denote the correlations between the visibilities at different baselines and frequencies. It should be noted that although we 
have shown S 2 and S3 as explicit functions of only the frequency differences Av, all these correlations also depend on the the 
central value v which is not shown as an explicit argument. Further, throughout our analysis we assume that all frequency 
differences are much smaller than the central frequency ie. Av/v <C 1. 

The correlation S 2 (Ui, U 2 , Av) between the visibilities at two baselines and frequencies has been calculated earlier 
(Bharadwaj & Sethi 2001; Bharadwaj & Pandey 2003; BA) who find that S 2 ~ if U 2 7^ — Ui. This is a consequence of the 
statistical homogeneity of the HI fluctuations. It is sufficient to restrict the analysis to Ui = — U 2 = U which we denote as 
S 2 (U, Av), and we have 

S 2 (U, Av) = f dfc„ P m (k) cos(k n r' v Au) . (12) 

lrv Jo 

were k = feym + (2n/r v )\J and r v = dr v /dv. The vector k has components fey and (2ir/r v )U respectively parallel and 
perpendicular to the line of sight. The fact that Pm(k), which includes redshift distortion, is isotropic in the directions 
perpendicular to the line of sight implies that S 2 is isotropic in U and we can write S 2 (U,Av). We also note that S 2 is real 
for the HI signal. This follows from the fact that Pm(k) is real and it is unchanged if fey — > —fey. 

The correlation of the visibilities at three different baselines and frequencies, S3 is the quantity of interest in this paper. 
This will be related to the HI bispectrum, Here, as for the power spectrum, we assume that Av/v <g; 1 , whereby the only 
term in eq. (9) for the visibility V(U, v + Av) which is affected by v — > v + Av is e - ,k \\ rv +* u , which can be approximated as 

g — ifcy (r v +r v Av) 

We then have 

S 3 (Ui, U 2 , U3, Aia, Av 2 ) = / d 3 fci d 3 k 2 d 3 k 3 o(Ui - ^k lx )a(U 2 - ^k 2± ) x 



(27r) 6 J ° y 2tt ^' v 2tt 

a(U 3 - ^k 3 ±) e - l ( fc ni+ fc 2|i+fe3||)^ e -i(*m A, 1+ fe 2|l A, 2 )r' v ^ (ki + k2 + k;j) Sm(kl)k2;k3) (13) 

It is convenient to write the d 3 k integrals as dk\\ d 2 k± and integrate over dfc 3 y, whereby the term e _l ( fc i||+ fe 2||+ fe 3||) r i' drops out 
because of the Dirac delta function. Also, we introduce a new variable y = k — (27r/r„)U and use the explicit form for the 
function o(U), whereby we have 

S 3 (Ul,U 2 ,U3,Al/l,Al*) = J^yr J dk n dk n e -i(fclH A. 1+ fc 2| | A, 2 ).: J ^ ^ ^ ^ 

*2,[(27r/r„)(Ui +U 2 + U 3 ) + yi + y 2 + y 3 ](7r6>§) 3 exp[-( r ^ /2) 2 (^ + yl + yl))B m (14) 

where the arguments of the bispectrum change as we carry out the integrals, but we do not show them explicitly. 
Carrying out the d 2 yz integral we have 

S 3 (Ui 1 U a ,U s> Avi,A^) = ^5 J dk n dk n e -i(fcm A. 1+ fc 2| | a, 2) .: ^S yi d 2 y 2 x 

exp[-(rvft,/2) 2 (y? + y\)\ cxp[-(r^ /2) 2 {yi + y 2 + (27r/r„)(Ui + U 2 + U 3 )} 2 ] Shi (15) 

The point to note is the two Gaussian functions exp[— (r v 8o/2) 2 (y( + y 2 )] and exp [-(r v e /2) 2 { yi +y 2 + (2tt/^)(Ui + U 2 + 
U3)} 2 ] are peaked around different values of yi and y 2 . While the former is peaked around yi = y 2 = 0, the latter is peaked 
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around yi + y2 = (— 27r/r„)(Ui + U2 + U3). The peaks of the two functions have very little overlap if | Ui + U2 + U3 |> 0, 
and the visibility correlations are exponentially suppressed if the vector sum of the baselines differs from zero. There are 
substantial correlations only for the sets of baselines for which | Ui + U2 + U3 < (irOo) 1 - In the rest of our analysis we 
only consider combinations of baselines for which Ui + U2 + U3 = 0, and the product of the two Gaussian functions becomes 
exp[— 2 (r„#o/2) 2 (j/i + y\ + yi ■ y2)]. This can be further simplified if the baselines we are dealing with are much larger 
than 1/ (7r0o )■ We can then approximate this function by a product of two Dirac delta functions « (167r/3)(r„#o) -4 5f)(yi + 
y2/2)<Jfj(y2). Using this in eq. (15) we have 

S 3 (Ui,U 2 ,U 3 ,AzA,Az, 2 ) = -|^r f dk n dk n e"^ 1 " A " 1+fc »H A " 2) ^ B m (ki, k 2 , k 3 ) (16) 

12 7r rS J 

where ki = k^m + (27r/r„)Ui, k2 = fc 2 ||m + (2-7r/r„)U2 and k ;j = — (fc^ + fc 2 ||)m + (2-7r/r„)U3. Further, it can be verified 
that S3 is real, and 

j3 q2 r ( 

S 3 (Ui,U2,U 3 ,Au 1 ,Av 2 ) = " ° 4 / dfci|| dk n cospiy Av x + k 2 \\ Av 2 )r„] Bm(ki, k 2 , k 3 ) (17) 

where we have also incorporated the fact that S3 depends only the triangle formed by Ui,U 2 and U3 which is completely 
specified by just the magnitudes (Ui, U2, U3). 

We use eqs. (12) and (17) to calculate the visibility correlations expected during the epoch of reionization. 



3 A MODEL FOR THE HI DISTRIBUTION 

The reionization of the HI in the universe started, possibly at a redshift z ~ 30, when the first luminous objects were formed. 
The radiation from these luminous objects and from the subsequently formed luminous objects ionized the low density HI in 
the universe. The reionization commences in small spherical regions (Stromgren sphere) surrounding the luminous objects. 
These spheres are filled with ionized HII gas, the rest of the universe being filled with HI . Gradually these ionized regions 
grow until they finally overlap, filling up the whole of space, and all the low density gas in the universe is ionized. The HI 
distribution during reionization is largely determined by the ionized regions. This is expected to be highly non-Gaussian 
carrying signatures of the size, shape and distribution of the discrete ionized regions. Here we adopt a simple model for the 
ionized regions. Though simple, this model suffices to illustrate the non-Gaussian nature of the HI distribution and allows us 
to calculate some of the salient observable consequences. 

We assume that the HI gas is heated well before it is reionized, and that the spin temperature is coupled to the gas 
temperature with T a T 7 so that (1 — T 7 /T s ) — » 1. It then follows that rim > (eq. 2) ie. the HI will be seen in emission. 
At any epoch a fraction of the volume fy is completely ionized, the ionized gas being in non-overlapping spheres of comoving 
radius R, the centers of the spheres being randomly distributed. This model is similar to that used by ZFH in the context 
of HI emission, and Gruzinov & Hu (1998) and Knox at al. (1998) in the context of the effect of patchy reionization on the 
CMBR. One would expect the centers of the ionized spheres to be clustered, given the fact that we identify them with the 
locations of the first luminous objects which are believed to have formed at the peaks of the density fluctuations. This effect, 
included in BA, has not been taken into account here. 

Following ZFH, we assume that the mean neutral fraction xm at any epoch is given by 

xm(z) = — jj^- (18) 

1 +cxp((z - z )/Az) 

with zo — 10 and Az — 0.5 so that 50% os the hydrogen is neutral at a redshift z = 10. The mean comoving number density 
of ionized spheres nm is related to the quantities defined earlier as fy = 1 — xm = (4-7ri? 3 /3)nHi. We have kept R as a free 
parameter and have used this to determine ram- 

We assume that the total hydrogen density traces the dark matter and hence it is ph(1+<5) where S refers to the fluctuations 
in the dark matter distribution. Then, in our model, the HI density is Phi(x, 2) = Ph(1 + S) [l — ^ a 9(\ x — x a j /R)] , where 
a refers to the different ionized spheres with centers at x a , and 9(y) is the Heaviside step function defined such that 9(y) = l 
for < y < 1 and zero otherwise. We then have 



r?m(x,z) = 



1 + zdv 



H(z) dr 



(19) 



where v refers to the peculiar velocity caused by 5. The point to note is that 77m (x, z) has contributions from two distinct effects 
namely the fluctuations arising from the gravitational clustering of the hydrogen which follows the dark matter distribution 
and the discrete ionized regions. Earlier studies (ZFH) have shown that the contribution from the discrete ionized regions 
dominates the HI power spectrum on length-scales larger than the size of the individual ionized bubbles at redshifts z ~ 10 
when fy ~ 0.5 and the HI signal is expected to be maximum. In the standard scenario, the initial dark matter fluctuation 
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8 is e assumed to be a Gaussian random field for which the bispectrum is zero. Non-Gaussian features of order ~ 8 2 arise 
from non-linear effects as the density fluctuation grow, but these effects are expected to be very small on the length scales 
of our interest at redshifts z > 10. The bispectrum Shi too will be dominated by the non-Gaussian features arising from the 
discrete ionized regions. Further, we expect the gravitational clustering of the hydrogen to make a smaller contribution to 
the bispectrum than it does to the power spectrum. The aim here being to investigate the non-Gaussian effects through a 
study of the bispectrum, it is justified to focus on just the contribution arising from the individual ionized regions, ignoring 
the effects of gravitation clustering. Under this assumption 



r?m(x,z) = 



(20) 



and it Fourier transform for k > is 
-fv W{kR) 



r?Hi(k,2) = > e (21) 

tihi * — ' 

a 

where W(y) = (S/y s )[s'm(y) — ycos(y)] is the spherical top hat window function. Using these we have 

= tiwHkR) 

tihi 

and 

BM^M) = - f v W{klR)W 2 {k * R)W{k3R) (23) 

n HI 

respectively for the power spectrum and the bispectrum. We use these to calculate the visibility correlations expected in this 
model. 

Our model has a limitation that it cannot be used when a large fraction of the volume is ionized as the ionized spheres 
start to overlap and the HI density becomes negative in the overlapping regions. Calculating the fraction of the total volume 
where the HI density is negative, we find this to be f v /2. We use this to asses the range of validity of our model. We restrict 
the model to z > 10 where f v < 0.5, and the HI density is negative in less than 12.5% of the total volume. 



4 RESULTS AND DISCUSSION 

In this section we present results for the visibility correlations expected from HI during the epoch of reionization. Our aim 
being to illustrate the non-Gaussian nature of the expected signal and its dependence on the ionized regions, we show results 
centered on only at a single frequency namely 125 MHz. This corresponds to a redshift z — 10.4 when the mean neutral 
fraction is xm = 0.67 (ie. fv = 0.33). We choose this particular frequency as xm is quite close to 0.5 where the HI signal 
is expected to be maximum, simultaneously ensuring that the volume fraction where the HI density predicted by our model 
becomes negative is small (~ 5%). Further, the HI signal is expected to be dominated by discrete ionized regions and hence 
we anticipate significant non-Gaussian features. 

We have used eqs. (12) and (17) to calculate the expected correlations between two and three visibilities respectively. 
For this it is necessary to specify a value for do, the beam size of the individual antennas in the array. Further, it may be 
noted that 8o « 0.6 x #fwhm- The value of 8o will depend on the physical dimensions of the antennas and the wavelength 
of observation. For the GMRT 8 = 1° at 325 MHz. We scale this using 8 oc A to obtain 8 = 2.6° at 125 MHz which we 
use here. The HI signal predicted here is for observations using the GMRT, and they can be directly compared to those in 
BA. Both S2 and S3 scale as 8q, and it is straightforward to scale the results presented here to make visibility correlation 
predictions for other radio telescopes. 

The comoving radius of the ionized spheres R is a free parameter in our model. Investigations on the growth of the 
ionized spheres (FZH) show that these will be at least a few Mpc in radius (possibly larger) at the redshift of interest. 
We have considered three possible values R = (1, 3, 5) hT 1 Mpc for which the respective values of rim are (78, 2.9, 0.63) x 
10- 3 /i 3 Mp C - 3 . 

For ease of graphical presentation, we have restricted our analysis of S3 to equilateral triangles for which the size of the 
baseline U completely specifies the triangle, and we have S3 ((7, Av). Further, we first consider the correlations at the same 
frequency ie. Av = Figures (1) and (2) show the results [S 2 ((7)] 1/2 and [-S 3 (t/)] 1/3 respectively. 

We find that at small U, [S2(U)] 1 ^ 2 is more or less constant with a value of the order of ~ 0.2 mjy for R = Sh^ 1 Mpc. 
The signal is proportional to R 3 ^ 2 and its magnitude increases as the ionized spheres become larger. Each baseline U can be 
associated with a comoving length-scale r v /(2nU) at the comoving distance where the HI radiation originated. The signal 
from the ionized spheres is constant across the baselines for which r 1/ /(2nU) is larger than the size of the spheres, and 
the signal falls at baselines for which U > r v /(2nR). Each baselines resolves out features larger than r^/(2-wU), and the 
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Figure 1. This shows the visibility correlation S2(U, Av) 1 / 2 as a function of U for Av = 0, for different values of R, the comoving radius 
of the ionized spheres. These predictions are for observations centered at 125 MHz. 
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Figure 2. This shows the visibility correlation [— S3 ((7, Au)] 1 ^ 3 as a function of U for Av = 0, for different values of R, the comoving 
radius of the ionized spheres. These predictions are for observations centered at 125 MHz. 

presence of discrete ionized regions make very little contribution to the signal at the large baselines. Comparing the results 
for [S2(£/)] 1//2 presented here with those presented in BA which also includes the effects of gravitational clustering, we note 
that the gravitational clustering signal is also of the order ~ 0.1 mjy at small U. The gravitational clustering signal also falls 
with increasing U , and the combined signal would depend critically on the size of the bubbles. For example, the signal from 
discrete ionized sources would dominate over the gravitational clustering signal at baselines U > 500 if the ionized spheres 
had comoving radius R = 5/i~ Mpc, whereas the gravitational clustering signal would possibly dominate throughout for 
R= l/i^Mpc. 

Turning our attention next to S3 (Figure 2), the first point to note is that this is negative. The shape of S3 as a function 
of U is very similar to that of S2, and its magnitude is around [— S3((7)] 1//3 ~ 0.06 mjy at R = Mpc, which is around 

three times smaller than [S2(?7)] 1 ^ 2 . At small U, [— Ss(?7)] 1 ^ 3 is more or less constant. Although our results are restricted to 
equilateral triangles, we expect the correlations to be nearly constant for triangles of all shapes provided all the baselines 
satisfy U < r l/ /(2nR). The signal is proportional to R 2 and its magnitude increases a little faster than that of [S2(U)] 1 > 12 
as R is increased. We expect the dark matter density fluctuations at z > 10 to be well in the linear regime on comoving 
length-scales ~ 10 hT 1 Mpc or larger, and the contribution to S3 from non-linear gravitational clustering is expected to be 
very small on these scales. It may be noted that the contribution to S3 from linear gravitational clustering is exactly zero in 
the standard scenario where the initial density fluctuations are a Gaussian random field. Further, the comoving length-scale 
10 ft. -1 Mpc corresponds to the baseline U ~ 100 and we expect the contribution from individual ionized spheres considered 
here to be the dominant signal at these baselines. 

We next consider the correlations between the visibilities at different frequencies. Again, for ease of graphical presentation 
we have restricted our analysis of S3 (17, Avi, AV2) to equilateral triangles with the added restriction that Avi = Av2 = Av, 
so we have Ss(U,Av). We have shown results only for R = 3 ft -1 Mpc, but a similar behaviour is expected for other values 
also. We find that both S2(U, Av) and S3(Z7, Av) fall rapidly, in nearly the same fashion independent of U, and are very close 
to zero by Av w 0.5 MHz. 
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Figure 3. This shows the visibility correlation S2(U, Av) 1 / 12 as a function of Av for the three different values of U shown in the figure. 
The comoving radius of the ionized spheres is assumed to be R = 3 /i -1 Mpc and the predictions are for observations centered at 125 MHz. 
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This shows the visibility correlation [— S3 (U, Av)] 1 / 3 as a function of Av for the three different values of U shown in the 
comoving radius of the ionized spheres is assumed to be R = 3/i _1 Mpc and the predictions are for observations centered at 



One of the main challenges in observing cosmological HI is to extract it from various contaminants which are expected 
to swamp this signal. The contaminants include Galactic synchrotron emission, free-free emission from ionizing halos (Oh 
& Mack 2003), faint radio loud quasars (DiMatteo. et al. 2002) and synchrotron emission from low redshift galaxy clusters 
(DiMatteo. et al. 2004) . Fortunately, all of these foregrounds have smooth continuum spectra and we expect their contribution 
to the visibilities to be correlated over large Av, whereas the HI contribution is uncorrelated beyond 1 MHz or less. It is, in 
principle, straightforward to fit the visibility correlations S2 and S3 at large Av and remove any slowly varying component 
thereby separating the contaminants from the HI signal. We also use this opportunity to note that this is a major advantage 
of using visibility correlations as compared to the angular power spectrum which exhibits substantial correlations even at two 
frequencies separated by ~ 10 MHz (Santos, Cooray & Knox 2004). 

An important fact which emerges from our analysis is that the HI signal in the correlation between three visibilities is 
of the same order as the correlation between two visibilities, the former being around three times smaller. This is a generic 
feature of the EOR HI signal, valid if the ionized regions are bubbles of the size R — lh^ 1 Mpc or larger. This signal arises 
from the Poisson noise of the discrete ionized regions, and it is enhanced if the size of the bubbles is increased. The fact that 
there is a substantial S3 tells us that there are large phase correlations between the visibilities. This is a consequence of the 
fact that there are only a few coherent features (the ionized regions) which dominate the whole HI signal. 

Investigations on the growth of ionized regions (FZH) show that there will be a spread in the sizes of the ionized regions 
at any given epoch. This will smoothen some of the sharp features seen in Figures (1) and (2). The ringing seen in the these 
figures is an artifact of there being only a single value of R and we do not expect this feature to be there if we have a spread in 
R. Further, the gravitational clustering signal not shown here may also dominate at large U. Despite all these limitations, we 
can still expect substantial correlations between three visibilities in a a more realistic analysis, this being a robust signature of 
the fact that reionization occurs through a few, large (R ~ a few Mpc) bubbles of ionized gas and the HI signal is dominated 
by Poisson noise. 

We next briefly discuss the noise levels and the integration times required to observe the HI signal, particularly addressing 
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the question whether S3 can be detected with integration times comparable to those needed for S2. We consider an array 
of N antennas, the observations lasting a time duration t, with frequency channels of width 8u spanning a total bandwidth 
B. It should be noted that the effect of a finite channel width 5v has not been included in our calculation which assumes 
infinite frequency resolution. This effect can be easily included by convolving our results for the visibility correlation with the 
frequency response function of a single channel. Preferably, 5v should be much smaller than the frequency separation at which 
the visibility correlation become uncorrelated. We use S to denote the frequency separation within which the visibilities are 
correlated, and beyond which they become uncorrelated. 

We use N2 and N3 to denote the rms. noise in S2 and S3 respectively. It is well known that N2 = ( ^ kB ^SYs ^ _i_ 

(Thompson, Moran & Swenson 1986), and we have N3 ~ ( 2k BTgvs ^ _J_^- assuming that we have Gaussian random 
noise, where Tsys is the system temperature and A e f is the effective area of a single antenna. The noise contributions will 
be reduced by a factor N if we combine N independent samples of the visibility correlation. A possible observational 
strategy for a preliminary detection of the HI signal would be to combine the visibility correlations at all baselines and 
frequency separations where there is a reasonable amount of signal. This gives N = [N(N — l)/2] (B/Sv) {S/8v) for the 
two visibility correlation and N = [N(N — l)(N — 2)/6] (B/Su) (S/8v) 2 for the three visibility correlations. It should be 
noted that we have used the fact that the S3 is non-zero only for the baselines between triplets of antennas. Combining all 

of this we have [N2] 1 / 2 ~ (»^) [ N{N % BS ] ^ ^ and [N3] 1 / 3 ~ (»^) [ N(N _ 1)( 6 N _ 2)BS , ] ^ ^y- The ratio 

[N 3 ] 1/3 /[N 2 ] 1/2 ~ [N(N - 1)B/(N - 2) 2 S'] 1/12 has a very weak dependence on N, B and S for a reasonable choice of values, 
and is of order unity. We thus see that, for a given integration time, we will achieve comparable noise levels in both the two 
and three visibility correlations. Estimates of the integration time to detect S2 (or equivalently the angular power spectrum) 
(BA, ZFH) indicate this to be around a few hundred hours. We find that it should be possible to also detect S3 in a comparable 
integration time. 
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